{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $p$ 次多项式空间的编号和指标的对应关系"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{array}{rccccccccccccc}\n",
    "p=0:&    &    &    &    &    &    &  1\\\\\n",
    "p=1:&    &    &    &    &    &  x &    &  y\\\\\n",
    "p=2:&    &    &    &    &  x^2 &    &  xy &    &  y^2\\\\\n",
    "p=3:&    &    &    &  x^3 &    &  x^2y &    &  xy^2 &    &  y^3\\\\\n",
    "p=4:&    &    &  x^4 &    &  x^3y &    &  x^2y^2 &    &  xy^3 &    &  y^4\\\\\n",
    "p=5:&    &  x^5 &    &  x^4y &    & x^3y^2 &    & x^2y^3 &    &  xy^4 &    &  y^5\\\\\n",
    "p=6:&  x^6 &    &  x^5y &    & x^4y^2 &    & x^3y^3 &    & x^2y^4 &    &  xy^5 &    &  x^6\\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二维情形"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "二重指标　$\\alpha=(a, b)$ , $\\mathbf x^\\alpha = x^{a}y^{b}$. $m_{\\alpha}$ 是 $\\alpha$ 所对应的单项式, $\\nabla{m_{\\alpha}}$ 是单项式的梯度, $\\Delta{m_{\\alpha}}$ 是单项式的拉普拉斯."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\alpha$ 所对应的编号为:\n",
    "\n",
    "$$\n",
    "0\\\\\n",
    "1 \\qquad 2\\\\\n",
    "3\\qquad 4 \\qquad 5\\\\\n",
    "6\\qquad 7 \\qquad 8\\qquad 9\\\\\n",
    "10\\qquad 11 \\qquad 12\\qquad 13\\qquad 14\\\\\n",
    "15\\qquad 16 \\qquad 17\\qquad 18\\qquad 19\\qquad 20\\\\\n",
    "\\cdots\\qquad \\cdots\\qquad \\cdots\\qquad \\cdots\\qquad \\cdots\\qquad \\cdots\\qquad\\cdots\\qquad\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\alpha$ 所对应的二重指标为:\n",
    "\n",
    "$$\n",
    "(0,0)\\\\\n",
    "(1,0) \\qquad (0,1)\\\\\n",
    "(2,0)\\qquad (1,1) \\qquad (0,2)\\\\\n",
    "(3,0)\\qquad (2,1) \\qquad (1,2)\\qquad (0,3)\\\\\n",
    "(4,0)\\qquad (3,1) \\qquad (2,2)\\qquad (1,3)\\qquad (0,4)\\\\\n",
    "(5,0)\\qquad (4,1) \\qquad (3,2)\\qquad (2,3)\\qquad (1,4)\\qquad (0,5)\\\\\n",
    "\\cdots\\qquad \\cdots\\qquad \\cdots\\qquad \\cdots\\qquad \\cdots\\qquad \\cdots\\qquad\\cdots\\qquad\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "若已知 $(a,b)$, 则对应的 $\\alpha$ 编号为:\n",
    "\n",
    "$$\n",
    "\\alpha = \\frac{(a+b)(a+b+1)}{2}+b\n",
    "$$\n",
    "\n",
    "若已知 $\\alpha$, 则对应的 $(a, b)$, 计算公式为\n",
    "$$\n",
    "a+b = \\lfloor \\frac{-1+\\sqrt{1+8\\alpha}}{2}\\rfloor\\\\\n",
    "b = \\alpha - \\frac{(a+b)(a+b+1)}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 三维情形"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "三重指标　$\\alpha=(a, b, c)$ , $\\mathbf x^\\alpha = x^{a}y^{b}z^{c}$. $m_{\\alpha}$ 是 $\\alpha$ 所对应的单项式, $\\nabla{m_{\\alpha}}$ 是单项式的梯度, $\\Delta{m_{\\alpha}}$ 是单项式的拉普拉斯."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "若已知 $(a,b)$, 则对应的 $\\alpha$ 编号为:\n",
    "\n",
    "$$\n",
    "\\alpha = \\frac{(a+b+c)(a+b+c+1)(a+b+c+2)}{6}+\\frac{(b+c)(b+c+1)}{2}+c\n",
    "$$\n",
    "\n",
    "若已知 $\\alpha$, 则对应的 $(a, b, c)$, 计算公式为\n",
    "$$\n",
    "a+b+c = \\left\\lfloor \\left(3\\alpha+\\frac{1}{3}\\sqrt{81\\alpha^2-\\frac{1}{3}}\\right)^{\\frac{1}{3}}+\\frac{1}{3\\left(3\\alpha+\\frac{1}{3}\\sqrt{81\\alpha^2-\\frac{1}{3}}\\right)^{\\frac{1}{3}}}-1\\right\\rfloor\\\\\n",
    "\\frac{(b+c)(b+c+1)}{2}+c = \\alpha - \\frac{(a+b+c)(a+b+c+1)(a+b+c+2)}{6}\n",
    "$$\n",
    "\n",
    "设 $\\alpha - \\frac{(a+b+c)(a+b+c+1)(a+b+c+2)}{6} = k$ ,则\n",
    "\n",
    "$$\n",
    "b+c = \\lfloor \\frac{-1+\\sqrt{1+8k}}{2}\\rfloor\\\\\n",
    "c = k - \\frac{(b+c)(b+c+1)}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 三维编号的转换关系"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "0\\\\\n",
    "1\\qquad 2\\\\\n",
    "3\\qquad 4 \\qquad 5\\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "四面体网格的每个单元的面上局部编号和全局编号之间的关系是什么呢？ 全局编号时每个面上从第一个顶点开始按三角形的方式编号， 但每个单元内局部编号的起点可能不一致，这样需要对这个面上的全局编号重新排序。 如上图中本来的排序为：\n",
    "\n",
    "给定四面体的一个面三个点的排序， \n",
    "$$\n",
    "(i, j, k)\n",
    "$$\n",
    "按和这个排序进行全局自由度的编号， 现给定另一个排序\n",
    "\n",
    "$$\n",
    "(i_1, j_1, k_1)\n",
    "$$\n",
    "\n",
    "它对应一个局部自由度的排序， 其中 $(i_1, j_1, k_1)$ 的排序方向和 $(i, j, k)$ 一定不一致吗？\n",
    "\n",
    "首先， 每个面的是从它左边单元中取出来的， 其法线方向指向左边单元的外部， \n",
    "\n",
    "    localFace = np.array([[1, 2, 3], [0, 3, 2]， [0, 1, 3], [0, 2, 1]])\n",
    "\n",
    "根据每个单元中的自由度编号顺序， 它的四个面的取法：\n",
    "\n",
    "    localFace = np.array([[1, 2, 3], [0, 2, 3], [0, 1, 3], [0, 1, 2]])\n",
    "    \n",
    "注意第 1 个面和第 3 个面的法向指向单元内部。\n",
    "\n",
    "给定一个单元， 取它对应的局部编号面 $(i_1, j_1, k_1)$， 及对应的全局面 $(i, j, k)$, 分两种情况讨论：\n",
    "\n",
    "1. 它是全局面对应的左边单元， 那么它的局部编号面要么和全局面完全一致， 要么排序相反\n",
    "2. 它是全局面对应的右边单元， 则局部第 0 个面和第 2 个面一定和全局面编号相反；第 1 个和第 3 个面的排序就可能是全局面排序一个同顺序的排列\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n",
      "0 \n",
      " [[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n",
      "[ 0.  1.  2.]\n",
      "1 \n",
      " [[1 0 0]\n",
      " [0 0 1]\n",
      " [0 1 0]]\n",
      "[ 0.  2.  1.]\n",
      "2 \n",
      " [[0 1 0]\n",
      " [1 0 0]\n",
      " [0 0 1]]\n",
      "[ 1.  0.  2.]\n",
      "3 \n",
      " [[0 0 1]\n",
      " [1 0 0]\n",
      " [0 1 0]]\n",
      "[ 2.  0.  1.]\n",
      "4 \n",
      " [[0 0 1]\n",
      " [0 1 0]\n",
      " [1 0 0]]\n",
      "[ 2.  1.  0.]\n",
      "5 \n",
      " [[0 1 0]\n",
      " [0 0 1]\n",
      " [1 0 0]]\n",
      "[ 1.  2.  0.]\n",
      "[[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n",
      "$$ \n",
      "\n",
      "0 \\\\ \n",
      "\n",
      "1 \\qquad  2 \\\\ \n",
      "\n",
      "$$\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def print_tri(p):\n",
    "    print('$$','\\n')\n",
    "    idx = -1\n",
    "    for i in range(p+1):\n",
    "        idx += 1\n",
    "        print(idx, end=' ')\n",
    "        for j in range(1, i+1):\n",
    "            idx += 1\n",
    "            print('\\\\qquad ', idx, end=' ')\n",
    "        print('\\\\\\\\', '\\n')\n",
    "    print('$$')\n",
    "        \n",
    "    \n",
    "def cell_idx(p):\n",
    "    ldof = int((p+1)*(p+2)/2)\n",
    "    idx = np.arange(0, ldof)\n",
    "    idx0 = np.floor((-1 + np.sqrt(1 + 8*idx))/2)\n",
    "    idx1 = idx0*(idx0+1)/2\n",
    "    cellIdx = np.zeros((ldof, 3), dtype=np.int)\n",
    "    cellIdx[:,2] = idx - idx1 \n",
    "    cellIdx[:,1] = idx0 - cellIdx[:,2]\n",
    "    cellIdx[:,0] = p - cellIdx[:, 1] - cellIdx[:, 2] \n",
    "    return cellIdx\n",
    "p = 1   \n",
    "idx = cell_idx(p)\n",
    "print(idx)\n",
    "localIdx = np.array([\n",
    "        [0, 1, 2],\n",
    "        [0, 2, 1],\n",
    "        [1, 0, 2],\n",
    "        [1, 2, 0],\n",
    "        [2, 1, 0],\n",
    "        [2, 0, 1],\n",
    "    ])\n",
    "for i in range(6):\n",
    "    I = idx[:, localIdx[i]]\n",
    "    print(i,'\\n', I)\n",
    "    J = I[:, 1] + I[:, 2]\n",
    "    print(J*(J+1)/2 +I[:, 2])\n",
    "print(idx)\n",
    "print_tri(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 0 0 0]\n",
      " [1 1 0 0]\n",
      " [1 0 1 0]\n",
      " [1 0 0 1]\n",
      " [0 2 0 0]\n",
      " [0 1 1 0]\n",
      " [0 1 0 1]\n",
      " [0 0 2 0]\n",
      " [0 0 1 1]\n",
      " [0 0 0 2]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "p = 2\n",
    "ldof = int((p+1)*(p+2)*(p+3)/6)\n",
    "idx = np.arange(1, ldof)\n",
    "idx0 = (3*idx + np.sqrt(81*idx*idx - 1/3)/3)**(1/3) \n",
    "idx0 = np.floor(idx0 + 1/idx0/3 - 1 + 1e-4) # a+b+c\n",
    "idx1 = idx - idx0*(idx0 + 1)*(idx0 + 2)/6\n",
    "idx2 = np.floor((-1 + np.sqrt(1 + 8*idx1))/2) # b+c\n",
    "cellIdx = np.zeros((ldof, 4), dtype=np.int)\n",
    "cellIdx[1:, 3] = idx1 - idx2*(idx2 + 1)/2\n",
    "cellIdx[1:, 2] = idx2 - cellIdx[1:, 3]\n",
    "cellIdx[1:, 1] = idx0 - idx2\n",
    "cellIdx[:, 0] = p - np.sum(cellIdx[:, 1:], axis=1)\n",
    "print(cellIdx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定两个相同形状的二维数组 `a` 和 `b`，其中 `b` 的每一行是 `a` 的对应行的一个排列， 请构造矩阵 `I`， 满足  `b = a[:, I]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[4 7 3]\n",
      " [7 1 8]\n",
      " [9 3 2]]\n",
      "[[7 3 4]\n",
      " [1 7 8]\n",
      " [9 2 3]]\n",
      "[[1 2 0]\n",
      " [1 0 2]\n",
      " [0 2 1]]\n",
      "[[7 3 4]\n",
      " [1 7 8]\n",
      " [9 2 3]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[[3, 7, 4],\n",
       "        [7, 3, 4],\n",
       "        [4, 7, 3]],\n",
       "\n",
       "       [[8, 1, 7],\n",
       "        [1, 8, 7],\n",
       "        [7, 1, 8]],\n",
       "\n",
       "       [[2, 3, 9],\n",
       "        [3, 2, 9],\n",
       "        [9, 3, 2]]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "a = np.array([\n",
    "        [4, 7, 3],\n",
    "        [7, 1, 8],\n",
    "        [9, 3, 2],\n",
    "    ])\n",
    "b = np.array([\n",
    "        [7, 3, 4],\n",
    "        [1, 7, 8],\n",
    "        [9, 2, 3],\n",
    "    ])\n",
    "print(a)\n",
    "print(b)\n",
    "idxb = np.argsort(b, axis=1)\n",
    "idxb1 = np.argsort(idxb, axis=1)\n",
    "idxa = np.argsort(a, axis=1)\n",
    "idx = np.array([0, 1, 2]).reshape(-1, 1)\n",
    "print(idxa[idx, idxb1])\n",
    "print(a[idx, idxa[idx, idxb1]])\n",
    "a[:, [[2, 1, 0], [1, 2, 0], [0, 1, 2]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3 7 4]\n",
      " [1 8 7]\n",
      " [9 2 3]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([\n",
    "        [4, 7, 3],\n",
    "        [7, 1, 8],\n",
    "        [9, 3, 2],\n",
    "    ])\n",
    "\n",
    "I = np.array([\n",
    "        [2, 1, 0],\n",
    "        [1, 2, 0],\n",
    "        [0, 2, 1],\n",
    "    ])\n",
    "\n",
    "b = a[np.arange(3).reshape(-1, 1), I]\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
